!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the program is the main program to calculate the supersonic rotational
!  flow field use Method of Characteristics
!  Ref. to 'Gas Dynamics' by Maurice J. Zucrow & Joe D. Hoffmann
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-08-15
! revised  :
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! call  :
!   subroutine :
!       VariableDef 
!       GasCons
!       Thermo
!       InterP
!       DirectW
!       InverseW
!       Wall
!       SymmetryP
!       FreeP
!       Transonic
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!    
!   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program MOC_IRRO
   use VariableDef
   implicit none
!     Define Input Variable
!   real*8::
!     Define Temporary Variable
   real*8,allocatable::xbarma(:)
   integer::i,j,k,k1,k2,k3,k4
   integer::flagexit,OnePN,OneP,Nshock,eoc,OverlapN,Overlap,next,shockexist
   integer,allocatable::ffn(:),csymmetry(:)
   real*8::thetai,alptemp,ytemp,atemp,stepk,stepi
!     Define Output Variable
   real*8,allocatable::x(:,:),y(:,:),Vu(:,:),Vv(:,:),theta(:,:),V(:,:),p(:,:),T(:,:),&
         &rho(:,:),M(:,:)
   
!********************************** program start *****************************************!
!****************    Initiate variable     ***
! open an file and write header
   open(unit=110,file='temp.txt')
   write(110,*) 'MOC Data'
   write(110,*)
   close(110)
! initiate the variables
   Tt=3000.0
   Pt=7000000.0
   g=1.2
   Rg=320.0
   p0=10000.0
   Vuu=1041.6
   pu=100000.0
   rhou=1.1613
   Tu=300.0
   Mu=3.0
   
   pi=3.141592653
   el=1e-5
   ep=1e-5
   erho=1e-5
   eV=1e-5
   etheta=1e-5
   delta=1
   icor=6
   
   yt=1.0 !0.1 !0.025 !1.0
   rtu=2.0 !0.2 !0.050 !2.0
   rtd=0.5 !0.05
   thetaa=15*pi/180
   thetae=10*pi/180
   xe=10.0 !1.0
   
   flagy40=0
   flagshock=0
   flagexit=0
   OnePN=0
   OneP=0
   Nshock=0
   OverlapN=0
   Overlap=0
   next=0
   shockexist=0
   
   Ny=11
   Nt=16
   Ni=99
   Nj=99
! allocate memory for array
   allocate(xbarma(Ny))
   allocate(ffn(Ni))
   allocate(csymmetry(Ni))
   
   allocate(x(Ni,Nj))
   allocate(y(Ni,Nj))
   allocate(Vu(Ni,Nj))
   allocate(Vv(Ni,Nj))
   allocate(theta(Ni,Nj))
   allocate(V(Ni,Nj))
   allocate(p(Ni,Nj))
   allocate(T(Ni,Nj))
   allocate(rho(Ni,Nj))
   allocate(M(Ni,Nj))
   
!***********************  initial line ******
   call Transonic(xbarma,x(1:Ny,1),y(1:Ny,1),Vu(1:Ny,1),Vv(1:Ny,1),theta(1:Ny,1),V(1:Ny,1)&
         &,M(1:Ny,1),T(1:Ny,1),p(1:Ny,1),rho(1:Ny,1))
   ffn(1)=1
   csymmetry(1)=1
   
!**********************   initial region *********
   do k1=2,Ny
      do k2=2,2*k1-2
         x1=x(k1,k2-1)
         y1=y(k1,k2-1)
         V1=V(k1,k2-1)
         theta1=theta(k1,k2-1)
         p1=p(k1,k2-1)
         rho1=rho(k1,k2-1)
         x2=x(k1-1,k2-1)
         y2=y(k1-1,k2-1)
         V2=V(k1-1,k2-1)
         theta2=theta(k1-1,k2-1)
         p2=p(k1-1,k2-1)
         rho2=rho(k1-1,k2-1)
         
         call InnerP()
         if (flagy40==1) then
            flagy40=0
            exit
         end if
         
         x(k1,k2)=x4
         y(k1,k2)=y4
         Vu(k1,k2)=Vu4
         Vv(k1,k2)=Vv4
         theta(k1,k2)=theta4
         V(k1,k2)=V4
         M(k1,k2)=M4
         T(k1,k2)=T4
         p(k1,k2)=p4
         rho(k1,k2)=rho4
      end do
      x1=x(k1,k2-1)
      y1=y(k1,k2-1)
      V1=V(k1,k2-1)
      theta1=theta(k1,k2-1)
      p1=p(k1,k2-1)
      rho1=rho(k1,k2-1)
      x3=x(k1-1,k2-2)
      y3=y(k1-1,k2-2)
      V3=V(k1-1,k2-2)
      theta3=theta(k1-1,k2-2)
      p3=p(k1-1,k2-2)
      rho3=rho(k1-1,k2-2)
      
      call SymmetryP()
      
      x(k1,k2)=x4
      y(k1,k2)=y4
      Vu(k1,k2)=Vu4
      Vv(k1,k2)=Vv4
      theta(k1,k2)=theta4
      V(k1,k2)=V4
      M(k1,k2)=M4
      T(k1,k2)=T4
      p(k1,k2)=p4
      rho(k1,k2)=rho4
      
      ffn(k1)=k2
      csymmetry(k1)=1
   end do
! write the initial region data
   open(unit=110,file='temp.txt',access='append',status='old')
   write(110,*) 'Initial Region Points'
   write(110,*) 'csymmetry i j x y Vu Vv theta V M T p rou'
   do i=1,Ny
      do j=1,2*i-1
         write(110,*) csymmetry(i),i,j,x(i,j),y(i,j),Vu(i,j),Vv(i,j),theta(i,j),V(i,j),&
               &M(i,j),T(i,j),p(i,j),rho(i,j)
      end do
   end do
   write(110,*)
   write(110,*) 'Mass flow rate and Thrust for initial line'
   write(110,*) 'mdot F mdot1d Fs F1d CDm etaF etaI Isp Isp1d'
   write(110,*) mdot,F,mdot1d,Fs,F1d,CDm,etaF,etaI,Isp,Isp1d
   write(110,*)
   write(110,*) 'Inverse Wall Region Points'
   write(110,*) 'csymmetry i j x y Vu Vv theta V M T p rou'
   close(110)
   
!****************************** inverse wall region   ************
   thetai=thetaa/(Nt-1)   ! initiate the step angle
   thetaa=0.0   ! initiate the wall angle
   k=2*Ny-1   ! initiate the first inverse wall c-line point number
   do k1=Ny+1,Ny+Nt-1
      thetaa=thetaa+thetai	!get the wall angle used in the cycle
      ! there the next c+ characteristics line must not intersect the wall before the 
      ! prescribed point
      do k3=2,k
         x(k1,1)=rtd*sin(thetaa)
         y(k1,1)=yt+rtd*(1-cos(thetaa))
         theta(k1,1)=thetaa
         
         x3=x(k1-1,1)
         y3=y(k1-1,1)
         V3=V(k1-1,1)
         theta3=theta(k1-1,1)
         p3=p(k1-1,1)
         rho3=rho(k1-1,1)
         x4=x(k1,1)
         y4=y(k1,1)
         theta4=theta(k1,1)
         
         alptemp=atan(Vv(k1-1,k3)/Vu(k1-1,k3))+asin(1/M(k1-1,k3))
         ytemp=y(k1-1,k3)+(x(k1,1)-x(k1-1,k3))*tan(alptemp)
         
         if (ytemp < y(k1,1)) then
            x1=x(k1-1,k3)
            y1=y(k1-1,k3)
            V1=V(k1-1,k3)
            theta1=theta(k1-1,k3)
            p1=p(k1-1,k3)
            rho1=rho(k1-1,k3)
            
            call IWall()
            
            Vu(k1,1)=Vu4
            Vv(k1,1)=Vv4
            V(k1,1)=V4
            M(k1,1)=M4
            T(k1,1)=T4
            p(k1,1)=p4
            rho(k1,1)=rho4
            exit
         end if
      end do
      
!      calculate thrust act on the wall point
      call Thrust(p(k1,1),y(k1,1),M(k1,1),p(k1-1,1),y(k1-1,1))
    
!     successive points
      do k2=2,k-k3+2
         x1=x(k1,k2-1)
         y1=y(k1,k2-1)
         V1=V(k1,k2-1)
         theta1=theta(k1,k2-1)
         p1=p(k1,k2-1)
         rho1=rho(k1,k2-1)
         x2=x(k1-1,k2+k3-2)
         y2=y(k1-1,k2+k3-2)
         V2=V(k1-1,k2+k3-2)
         theta2=theta(k1-1,k2+k3-2)
         p2=p(k1-1,k2+k3-2)
         rho2=rho(k1-1,k2+k3-2)
         
         call InnerP()
         if (flagy40==1) then
            flagy40=k2
            exit
         end if
         
         x(k1,k2)=x4
         y(k1,k2)=y4
         Vu(k1,k2)=Vu4
         Vv(k1,k2)=Vv4
         theta(k1,k2)=theta4
         V(k1,k2)=V4
         M(k1,k2)=M4
         T(k1,k2)=T4
         p(k1,k2)=p4
         rho(k1,k2)=rho4
      end do
      k=k2
!      axial point
      x1=x(k1,k2-1)
      y1=y(k1,k2-1)
      V1=V(k1,k2-1)
      theta1=theta(k1,k2-1)
      p1=p(k1,k2-1)
      rho1=rho(k1,k2-1)
      x3=x(k1-1,ffn(k1-1))
      y3=y(k1-1,ffn(k1-1))
      V3=V(k1-1,ffn(k1-1))
      theta3=theta(k1-1,ffn(k1-1))
      p3=p(k1-1,ffn(k1-1))
      rho3=rho(k1-1,ffn(k1-1))
      
      call SymmetryP()
      
      x(k1,k2)=x4
      y(k1,k2)=y4
      Vu(k1,k2)=Vu4
      Vv(k1,k2)=Vv4
      theta(k1,k2)=theta4
      V(k1,k2)=V4
      M(k1,k2)=M4
      T(k1,k2)=T4
      p(k1,k2)=p4
      rho(k1,k2)=rho4
      
      ! add a point between last and inverse second point on Ny+1 r-c-line
      ! or the r-c-line where flagy40==0 before the last point
      if ((flagy40==ffn(k1-1)) .or. (k1==Ny+1)) then   ! add point, begin
         flagy40=0
         k=1
         ! translate 
         x(k1,k2+k)=x(k1,k2)
         y(k1,k2+k)=y(k1,k2)
         V(k1,k2+k)=V(k1,k2)
         theta(k1,k2+k)=theta(k1,k2)
         p(k1,k2+k)=p(k1,k2)
         rho(k1,k2+k)=rho(k1,k2)
         Vv(k1,k2+k)=V(k1,k2+k)*sin(theta(k1,k2+k))
         Vu(k1,k2+k)=V(k1,k2+k)*cos(theta(k1,k2+k))
         call Thermo(p(k1,k2+k),rho(k1,k2+k),V(k1,k2+k),atemp,T(k1,k2+k),M(k1,k2+k))
         ! interpolation
         stepk=0.0
         stepi=1.0/(k+1.0)
         do i=1,k
            stepk=stepk+stepi
            x(k1,k2+i-1)=x(k1,k2-1)+stepk*(x(k1,k2+k)-x(k1,k2-1))
            y(k1,k2+i-1)=y(k1,k2-1)+stepk*(y(k1,k2+k)-y(k1,k2-1))
            V(k1,k2+i-1)=V(k1,k2-1)+stepk*(V(k1,k2+k)-V(k1,k2-1))
            theta(k1,k2+i-1)=theta(k1,k2-1)+stepk*(theta(k1,k2+k)-theta(k1,k2-1))
            p(k1,k2+i-1)=p(k1,k2-1)+stepk*(p(k1,k2+k)-p(k1,k2-1))
            rho(k1,k2+i-1)=rho(k1,k2-1)+stepk*(rho(k1,k2+k)-rho(k1,k2-1))
            Vv(k1,k2+i-1)=V(k1,k2+i-1)*sin(theta(k1,k2+i-1))
            Vu(k1,k2+i-1)=V(k1,k2+i-1)*cos(theta(k1,k2+i-1))
            call Thermo(p(k1,k2+i-1),rho(k1,k2+i-1),V(k1,k2+i-1),atemp,T(k1,k2+i-1),&
                  &M(k1,k2+i-1))
         end do
         k2=k2+k
         k=k2
      end if   ! add point, end
      
      csymmetry(k1)=1   ! the c-line arrived axis
      ffn(k1)=k2
   end do

! write the data in the file
   open(unit=110,file='temp.txt',access='append',status='old')
   do i=Ny+1,Ny+Nt-1
      do j=1,ffn(i)
         write(110,*) csymmetry(i),i,j,x(i,j),y(i,j),Vu(i,j),Vv(i,j),theta(i,j),V(i,j),M(i,j),&
               &T(i,j),p(i,j),rho(i,j)
      end do
   end do
   write(110,*)
!   write data in file in the order of thrust, 1-d thrust, F/F1d, I/I1d, specific impulse, 1-d Isp
   write(110,*) 'thrust before direct wall'
   write(110,*) 'F F1d etaF etaI Isp Isp1d'
   write(110,*) F,F1d,etaF,etaI,Isp,Isp1d
   write(110,*)
   write(110,*) 'Direct Wall Region Points'
   write(110,*) 'csymmetry i j x y Vu Vv theta V M T p rou'
   close(110)
   
!***************************    direct wall region   ***************
   call Wall2c()     ! second order wall function
!   call Wall3c()   !third order wall function
   xa=x(k1-1,1)   ! xa,ya  attachment point location
   ya=y(k1-1,1)

   do k1=Ny+Nt,100
      x3=x(k1-OverlapN-1,1)
      y3=y(k1-OverlapN-1,1)
      V3=V(k1-OverlapN-1,1)
      theta3=theta(k1-OverlapN-1,1)
      p3=p(k1-OverlapN-1,1)
      rho3=rho(k1-OverlapN-1,1)
      if (OneP==1) then   ! if there the previous c-line contains only one point, and which on the wall
         OnePN=OnePN+1   ! the account that c-line contains only one point
         x2=x(k1-OverlapN-OnePN-1,OnePN+2+OverlapN)
         y2=y(k1-OverlapN-OnePN-1,OnePN+2+OverlapN)
         V2=V(k1-OverlapN-OnePN-1,OnePN+2+OverlapN)
         theta2=theta(k1-OverlapN-OnePN-1,OnePN+2+OverlapN)
         p2=p(k1-OverlapN-OnePN-1,OnePN+2+OverlapN)
         rho2=rho(k1-OverlapN-OnePN-1,OnePN+2+OverlapN)
      else
         OnePN=0
         x2=x(k1-OverlapN-1,2+OverlapN)
         y2=y(k1-OverlapN-1,2+OverlapN)
         V2=V(k1-OverlapN-1,2+OverlapN)
         theta2=theta(k1-OverlapN-1,2+OverlapN)
         p2=p(k1-OverlapN-1,2+OverlapN)
         rho2=rho(k1-OverlapN-1,2+OverlapN)
      end if
      OneP=0
      
      call DWall()
      
!      justify if there a point located before the front point on the c-line, overlapped
      if (x4 < x2) then
         x(k1,1)=0.0   ! write in flag point for delete 
         y(k1,1)=0.0
         Overlap=1   ! overlap occurred
         OverlapN=OverlapN+1   ! overlapped point number        
         cycle
      end if
      
      Overlap=0   ! overlap does not occurred
      
      x(k1,1)=x4
      y(k1,1)=y4
      Vu(k1,1)=Vu4
      Vv(k1,1)=Vv4
      theta(k1,1)=theta4
      V(k1,1)=V4
      M(k1,1)=M4
      T(k1,1)=T4
      p(k1,1)=p4
      rho(k1,1)=rho4
      k3=3   ! if the next point calculated directly from the corresponding previous c-line point
!      judge if it exceed the lip point, and 
!      the exit lip point calculated use IWall subroutine
!      this section analogy to that in the inverse wall section
      if (x4 > xe) then
         flagexit=1
         do k3=2,k
            x(k1,1)=xe
            y(k1,1)=ye
            
            x3=x(k1-OverlapN-1,1)
            y3=y(k1-OverlapN-1,1)
            V3=V(k1-OverlapN-1,1)
            theta3=theta(k1-OverlapN-1,1)
            p3=p(k1-OverlapN-1,1)
            rho3=rho(k1-OverlapN-1,1)
            x4=x(k1,1)
            y4=y(k1,1)
            theta4=thetae
            
            alptemp=atan(Vv(k1-OverlapN-1,k3)/Vu(k1-OverlapN-1,k3))+asin(1/M(k1-OverlapN-1,k3))
            ytemp=y(k1-OverlapN-1,k3)+(x(k1,1)-x(k1-OverlapN-1,k3))*tan(alptemp)
         
            if (ytemp < y(k1,1)) then
               x1=x(k1-OverlapN-1,k3)
               y1=y(k1-OverlapN-1,k3)
               V1=V(k1-OverlapN-1,k3)
               theta1=theta(k1-OverlapN-1,k3)
               p1=p(k1-OverlapN-1,k3)
               rho1=rho(k1-OverlapN-1,k3)
               
               call IWall()
               
               Vu(k1,1)=Vu4
               Vv(k1,1)=Vv4
               theta(k1,1)=theta4
               V(k1,1)=V4
               M(k1,1)=M4
               T(k1,1)=T4
               p(k1,1)=p4
               rho(k1,1)=rho4
               k=k-k3+3
               exit
            end if
         end do 
      end if
!      calculated thrust data
      call Thrust(p(k1,1),y(k1,1),M(k1,1),p(k1-1,1),y(k1-1,1))
      
!      successive points
      if(flagshock==1) then   ! previous c-line crossed
         Nshock=Nshock+1   ! the account of c-line crossed
      else
         Nshock=0
      end if
      flagshock=0   ! not crossed
      eoc=1   ! c-line ended at the axis
      Overlap=0
      do k2=2,k-1
! if overlap occurred, this point must calculated from the point which
! does not overlapped the front point on the same c-line, next
         if (next==1) then
            x1=x(k1,k2-1-Overlap)
            y1=y(k1,k2-1-Overlap)
            V1=V(k1,k2-1-Overlap)
            theta1=theta(k1,k2-1-Overlap)
            p1=p(k1,k2-1-Overlap)
            rho1=rho(k1,k2-1-Overlap)
         else
            x1=x(k1,k2-1)
            y1=y(k1,k2-1)
            V1=V(k1,k2-1)
            theta1=theta(k1,k2-1)
            p1=p(k1,k2-1)
            rho1=rho(k1,k2-1)
         end if
         next=0
         
!         this point calculated from the point on the previous Nshoch-th c-line point 
         do k4=0,Nshock
            if (k2+k3-2+k4+OverlapN <= ffn(k1-OverlapN-k4-1)) then
!            if the point just located after the shock (assumed as a c-line), it must be calculated
!            from the parameters after the shock
               if ((Nshock>0) .and. (k4==Nshock) .and. (shockexist==1)) then
                  x2=x(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  y2=y(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  V2=V(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  theta2=theta(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  p2=p(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  rho2=rho(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  ! call Oshock()
                  exit
               else
                  x2=x(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  y2=y(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  V2=V(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  theta2=theta(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  p2=p(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  rho2=rho(k1-OverlapN-k4-1,k2+k3-2+k4+OverlapN)
                  exit
               end if
            end if
         end do
         
         call InnerP()
         
         if (x4 < x1) then   ! point overlapped
            x(k1,k2)=0.0
            y(k1,k2)=0.0
            if (k2==k-1) exit   ! if it is the point before symmetry point,exit without Overlap and next values
            Overlap=Overlap+1
            next=1
            cycle
         end if
         
         if (flagy40==1) then   ! y<0, exit
            flagy40=0
            k=k2
            exit
         end if         
!   if there is a shock occured, MOC disabled and the current characteristics line
!   crossed the spatial front characteristics line, this c-line must be terminated, 
!   and the next characteristics line proceeding
         if (flagshock==1) then
            if (k2==2) OneP=1   ! if this point is the second point on the c-line
            shockexist=1
            eoc=0   ! c-line does not arrived axis
            exit
         end if
!         for skip the points that get close to axis too much and weird
         if (y4 < yt*0.01) exit
!         store the results
         x(k1,k2)=x4
         y(k1,k2)=y4
         Vu(k1,k2)=Vu4
         Vv(k1,k2)=Vv4
         theta(k1,k2)=theta4
         V(k1,k2)=V4
         M(k1,k2)=M4
         T(k1,k2)=T4
         p(k1,k2)=p4
         rho(k1,k2)=rho4
      end do
      
!      c-line terminated on the inter region, initiate next c-line
      if (eoc==0) then
         ffn(k1)=k2-1   ! account of point(s) on the c-line
         k=ffn(k1-Nshock-1)-k3+2-Nshock-Overlap   ! next c-line maximum point number
         csymmetry(k1)=0   ! the c-line does not crrived axis
         
         do i=1,ffn(k1)   ! delete the point(s) overlapped the front point, x=y=0
            if (x(k1,i)==0.0) then
               do j=i+1,ffn(k1)
                  if (x(k1,j) > 0) then
                     x(k1,i)=x(k1,j)
                     y(k1,i)=y(k1,j)
                     Vu(k1,i)=Vu(k1,j)
                     Vv(k1,i)=Vv(k1,j)
                     theta(k1,i)=theta(k1,j)
                     V(k1,i)=V(k1,j)
                     T(k1,i)=T(k1,j)
                     p(k1,i)=p(k1,j)
                     rho(k1,i)=rho(k1,j)
                     x(k1,j)=0.0
                     exit
                  end if
               end do
            end if
         end do
         
         ffn(k1)=ffn(k1)-Overlap   ! account of point(s) on the c-line
         k=k-Overlap   ! next c-line maximum point number
         if (flagexit==1) exit	! if the exit lip point arrived
         cycle
      end if
      
      k=k2
!      axial point
      j=0
      do i=k1-1,1,-1   ! get previous axial point
         if (csymmetry(i)==0) then
            j=j+1
         else
            exit
         end if
      end do
      x3=x(k1-j-1,ffn(k1-j-1))
      y3=y(k1-j-1,ffn(k1-j-1))
      V3=V(k1-j-1,ffn(k1-j-1))
      theta3=theta(k1-j-1,ffn(k1-j-1))
      p3=p(k1-j-1,ffn(k1-j-1))
      rho3=rho(k1-j-1,ffn(k1-j-1))
      if (next==1) then
         x1=x(k1,k2-1-Overlap)
         y1=y(k1,k2-1-Overlap)
         V1=V(k1,k2-1-Overlap)
         theta1=theta(k1,k2-1-Overlap)
         p1=p(k1,k2-1-Overlap)
         rho1=rho(k1,k2-1-Overlap)
      else
         x1=x(k1,k2-1)
         y1=y(k1,k2-1)
         V1=V(k1,k2-1)
         theta1=theta(k1,k2-1)
         p1=p(k1,k2-1)
         rho1=rho(k1,k2-1)
      end if
      
      call SymmetryP()
      
!      c-line crossed just before the symmetry point
      if ((k2 >= 2) .and. (x4 <= x(k1-Nshock-1,k2+k3-2+Nshock))) then
         flagshock=1
         ffn(k1)=k2-1
         k=ffn(k1-Nshock-1)-k3+2-Nshock-Overlap
         csymmetry(k1)=0
         
         do i=1,ffn(k1)
            if (x(k1,i)==0.0) then
               do j=i+1,ffn(k1)
                  if (x(k1,j) > 0) then
                     x(k1,i)=x(k1,j)
                     y(k1,i)=y(k1,j)
                     Vu(k1,i)=Vu(k1,j)
                     Vv(k1,i)=Vv(k1,j)
                     theta(k1,i)=theta(k1,j)
                     V(k1,i)=V(k1,j)
                     T(k1,i)=T(k1,j)
                     p(k1,i)=p(k1,j)
                     rho(k1,i)=rho(k1,j)
                     x(k1,j)=0.0
                     exit
                  end if
               end do
            end if
         end do
         
         ffn(k1)=ffn(k1)-Overlap
         k=k-Overlap
         if (flagexit==1) exit
         cycle
      end if
      
!      store the results
      x(k1,k2)=x4
      y(k1,k2)=y4
      Vu(k1,k2)=Vu4
      Vv(k1,k2)=Vv4
      theta(k1,k2)=theta4
      V(k1,k2)=V4
      M(k1,k2)=M4
      T(k1,k2)=T4
      p(k1,k2)=p4
      rho(k1,k2)=rho4
!      record point number along the c-line and flag if the c-line arrived axis
      ffn(k1)=k2
      csymmetry(k1)=1
      shockexist=0

      do i=1,ffn(k1)
         if (x(k1,i)==0.0) then
            do j=i+1,ffn(k1)
               if (x(k1,j) > 0) then
                  x(k1,i)=x(k1,j)
                  y(k1,i)=y(k1,j)
                  Vu(k1,i)=Vu(k1,j)
                  Vv(k1,i)=Vv(k1,j)
                  theta(k1,i)=theta(k1,j)
                  V(k1,i)=V(k1,j)
                  T(k1,i)=T(k1,j)
                  p(k1,i)=p(k1,j)
                  rho(k1,i)=rho(k1,j)
                  x(k1,j)=0.0
                  exit
               end if
            end do
         end if
      end do
      
      ffn(k1)=ffn(k1)-Overlap
      k=k-Overlap
      if (flagexit==1) exit
   end do
   
!   write to file
   open(unit=110,file='temp.txt',access='append',status='old')
   do i=Ny+Nt,k1
      do j=1,ffn(i)
         write(110,*) csymmetry(i),i,j,x(i,j),y(i,j),Vu(i,j),Vv(i,j),theta(i,j),V(i,j),&
               &M(i,j),T(i,j),p(i,j),rho(i,j)
      end do
   end do
   write(110,*)
   write(110,*) 'thrust of nozzle'
   write(110,*) 'F F1d etaF etaI Isp Isp1d'
   write(110,*) F,F1d,etaF,etaI,Isp,Isp1d
   close(110)
   write(*,*) 'complete!'
   
end program MOC_IRRO
